
/* This file contains functions which implement a sequential algorithm to model the steady state equations for a metbolic network */

/*
// includes, system
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
*/

void mnSeq(int n, int p, float *vf, float *vr, float *vt, float *a, float *dm, float *m, float *kf, float *kr, float *S, float epsilon, float tolerance, float t) 
{

int i, j;
float fwdStoic[n * p];  //substrate stoichiometric values in forward direction - Sf
float revStoic[n * p];  //substrate stoichiometric values in reverse direction - Sr
unsigned int timeStep;

//forward and reverse stoichiometric matrices
for(i = 0; i < n * p; i++)
{
       if(S[i] < 0)
       {
           fwdStoic[i] = -1 * S[i];
           revStoic[i] = 0.0f;
       }
       else if(S[i] > 0)
       {
           fwdStoic[i] = 0.0f;
           revStoic[i] = S[i];  
       }
       else 
       { 
           fwdStoic[i] = 0.0f;
           revStoic[i] = 0.0f;  
       }
}


for(timeStep = 0; timeStep < 1000; timeStep++)
{
// loop through the fluxes and calculate the fwd and reverse flux vectors
    for(i = 0; i < n; i++)
    {
        vf[i] = kf[i];
        vr[i] = kr[i];
        for(j = 0; j < p; j++)
        { 
            vf[i] = vf[i] * powf(m[j],fwdStoic[(j * n) + i]); // indices are reversed in matlab code
            vr[i] = vr[i] * powf(m[j],revStoic[(j * n) + i]);
        }
    }

    // calculate the metabolic concentration change

    for(i = 0; i < p; i++)
    {
        dm[i] = 0;
        for(j = 0; j < n; j++)
        {
                dm[i] += S[(i * n) + j] * (vf[j] - vr[j]);
        }
    }

    // Compute new metabolite concentrations
    for(i = 0; i < p; i++)
    {
        m[i] = m[i] + (a[i] * dm[i] * t);
    }

    for(i = 0; i < p; i++)
    {
        // if concentration is negative, set it to a very small value.
        if(m[i] < 0)
        {
            m[i] = .00001;
        }
    }

    // Compute new volume scale factor
    for(i = 0; i < p; i++)
    {
        a[i] = epsilon + ( 1 / (abs(dm[i]) + epsilon));
    }
}



// Calculate total flux
for(i = 0; i < n; i++)
{
   vt[i] = vf[i] - vr[i];
}

}



